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The hot gas embedded in the large-scale structures in the Universe produces secondary fluctua- 
tions in the Cosmic Microwave Background (CMB). Because it is proportional to the gas pressure 
integrated along the line of sight, this effect, the thermal Sunyaev-Zel'dovich (SZ) effect, provides 
a direct measure of large-scale structure and of cosmological parameters. We study the statistical 
properties of this effect using both hydrodynamical simulations and analytical predictions from an 
extended halo model. The Adaptive Mesh Refinement scheme, used in the newly developed code 
RAMSES, provides a dynamic range of 4 order of magnitudes, and thus allows us to significantly 
improve upon earlier calculations. After accounting for the finite mass resolution and box size of the 
simulation, we find that the halo model agrees well with the simulations. We discuss and quantify 
O ■ the uncertainty in both methods, and thus derive an accurate prediction for the SZ power spectrum 

in the 10 2 < Z < 10 5 range of multipole. We show how this combined analytical and numerical 
approach is essential for accuracy, and useful for the understanding of the physical processes and 
scales which contribute to large-scale SZ anisotropies. 
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J> ■ I. INTRODUCTION 

\o ; 

The hot gas in the IGM induces distortions in the spectrum of the Cosmic Microwave Background (CMB) through 
the Sunyaev-Zel'dovich (SZ) effect |^,^]. This effect is proportional to the integrated electron pressure along the 
line-of-sight and can be directly measured by mapping the anisotropies of the CMB (see ||,|| for reviews). It acts 
as a foreground for the measurement of primary anisotropies of the CMB (see eg. |5| for a review), and is a direct 
probe of large-sale structure and of the distribution of baryons in the universe. This technique is now well established 
for the study of individual clusters of galaxies (eg. A more complete and unbiased statistical description of 

the SZ effect requires wide field surveys, and will soon be achieved by upcoming and future CMB experiments (see 
P,Bp3j and reference therein). In this perspective, accurate theoretical predictions of the statistics of the large-scale 
SZ effect are necessary. 

The statistics of SZ anisotropies have been studied using various methods. Scaramella et al. |0|, and more recently 
da Silva et al. flj|,[il| and Seljak, Burwell & Pen pQ] , have used numerical simulations to construct SZ maps and 
study their statistical properties. Persi et al. Jll| and Refregier et al. (l^] used instead a semi-analytical method, 
, ' consisting of computing the SZ angular power spectrum by projecting the 3-dimensional power spectrum of the gas 
pressure on the sky. Aghanim et al. |8), Atrio-Barandela & Miicket [[llj and Komatsu & Kitayama JlJ] used the Press 



& Schechter formalism 35 to predict the SZ power spectrum in various CDM models. In the same spirit, Cooray |21 j 
recently used a self-consistent extended halo model |3(],|2j| to make analytical predictions for the large-scale SZ effect. 
In a different approach, Zhang & Pen [^2| recently presented analytical predictions based on hierarchical clustering 
and non-linear perturbation theory. 

The accuracy of these different predictions are limited both by the limited dynamical range of the simulations and the 
lack of detailed comparison between the numerical and analytical methods. In this paper, we combine the approaches 
of Refregier et al. Jl9[ ] and of Cooray pl| to carefully compare predictions from numerical simulations to that from the 
extended halo model. We use a recently developed Adaptative Mesh Refinement (AMR) hydrodynamical code, called 
RAMSES (see Teyssier [|3j for a complete presentation), to compute the power spectrum of the gas pressure and of 
the Dark Matter (DM) density in a ACDM universe. The AMR simulations provide a dynamic range of 4 orders of 
magnitude and thus allow us to significantly improve upon earlier calculations. After ensuring the proper treatment 
of the DM power spectrum, the extended halo model allows us to predict the gas-pressure power spectrum without 
additional free parameter. We compare the predictions from both methods in detail, and pay particular attention to 
the limitations imposed by finite mass resolution and box size. We then compute the resulting SZ power spectrum 
and discuss the reliability and accuracy of our prediction. 

This paper is organized as follows. In § Q , we briefly describe the SZ effect and derive expressions for the integrated 



comptonization parameter and the SZ power spectrum. In §111, we describe the halo model for the DM and the 
baryons. In j pV| , we describe the numerical methods used for the RAMSES simulations. In §|y| we present the 
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comparative results for the DM power spectrum, the evolution of the density- weighted temperature, the pressure 
power spectrum and the SZ power spectrum. Our conclusions are summarized in §VL 



II. SUNYAEV-ZEL'DOVICH EFFECT 

The thermal SZ effect is produced from the inverse Compton scattering of CMB photons by the hot electrons in the 
IGM The first SZ observable is the mean comptonization parameter which can be measured from the distortion 

of the CMB Planck spectrum (see pq] for a review) . Following the conventions of |L9| , we find that it is given by 

i-Xi 

■ j d\T p a^ 2 . (1) 
Jo 

where x is the comoving distance, a is the expansion parameter and yo ~ 1.710 x 10~ 16 (jfj^j K _1 Mpc -1 , for a 

He mass fraction of 0.24 and assuming thermal equilibrium between the ions and electrons. The limiting distance Xi 
corresponds to the reionization redshift, for which recent theoretical studies suggest a value around Zi ~ 10 p^| . The 
mean density-weighted temperature of the gas is proportional to the average pressure and is defined as 

T p = ( P T)/{p), (2) 

where p and T are the density and temperature of the gas, respectively. 

The large-scale SZ effect can also be measured from the anisotropics that it induces on the CMB temperature. In the 
Rayleigh-Jeans (RJ) regime and in the small angle approximation, the angular power spectrum of these fluctuations 
is given by (see again |l^] for conventions) 

rXi „ /f \ 



Q~4^y d X T p P p ^- lX ja-' l T- z 1 (3) 

where r is the comoving angular diameter distance, and P p (k, x) is the 3-dimensional power spectrum of the pressure 
fluctuations 5p = (P — (P))/(P) at comoving distance x- I n the course of our work, we found that y is sensitive 
to the actual value of Zi, while the SZ power spectrum C; is almost insensitive to it. For definitiveness, we hitherto 
assume a fixed value for the reionization redshift of Zi = 10. 

III. HALO MODEL 
A. Dark Matter 

To model the dark matter as a collection of halos, we follow the approach of Seljak [^9| and of Ma & Fry [ p0| . 
Because its precise predictions depend on a number of assumptions, we here summarize the components of our halo 
model. The number of halos of mass M per unit comoving volumes and per unit mass at redshift z is given by the 
mass function 

dn p dv 

dM = MdM fi ^ (4) 
where p is the average matter density. The peak height v is defined as 

y{M,z) = -^— (5) 
<J(M, z) 

where cr(M, z) is the linear rms mass fluctuation in spheres of radius R given by M = AifpBc' /3 at redshift z. This 
quantity can be computed from the linear power spectrum Pn n (k, z), which we evaluate from the the BBKS transfer 
function |3l| (with the conventions of |3^]) evolved with the linear growth factor. The density threshold S c is equal 
to 1.68 in an Einstcin-De Sitter universe, and has a weak cosmology dependence which we compute using the results 
of Kitayama & Suto [Q. We adopt the standard Press-Schechter formalism |35|| , which dictates 

/(") = \jl e ^ 2/2 - (6) 
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In order for the halos to amount to the total mass density p, the mass function must be normalized as 

=pl dM dM M = l (7) 

This is formally satisfied by the Press & Schechter mass function, as soon as a(M, z) — > +00 for M — > 0. However, 
in CDM-like model this divergence at small M's is very slow, thus, practically, leaving a fraction of the mass in the 
background (i.e. not in collapsed halos). We account for this background mass by adding a constant to the mass 
function in our smallest mass bin (typically 10 6 A/ Q ). 

To each halo, we assign a Navarro, Frenk & White (NFW) mass density profile J3(|, which is given by 

p(r) = Ps u{r/r s ), u(x) = x' 1 (1 + x)~ 2 (8) 

where r is the comoving radius, r s and p s are the characteristic radius and density, respectively . These char- 
acteristic quantities can be written in terms of the virial radius r v within which the mean density contrast is 
8 V = 200. Identifying the mass M to the virial mass, we get M = 4Tfpd v r^/3 = 47r drr 2 p(r), yielding 

p s = ~pS v c 3 [ln(l + c) — c/(l + c)] 1 /3. In the previous expression, we have used the compactness parameter which is 
defined as 



c(M, z) =r v /r s . 

For the A CDM model we consider here, we adopt the functional form of c given by |32|,E!l| 



c(M,z)~g(z) 



M 



M*{z) 



-h(z) 



(9) 



(10) 



where g(z) ~ 10.3(1 + z) 3 and h(z) ~ 0.24(1 + z) °' 3 . Here, M*(z) is the non-linear mass scale defined by 
v{M*, z) = 1. The bias for the clustering of the halos is given by the Mo & White formalism BTj], namely by 



b(M, z) = l + 



v z - 1 



(11) 



Assuming that the matter distribution can be modeled as a collection of these halos, the matter power spectrum is 
given by p9p0| 



P(k,z) =P 1 (k,z) + P 2 (k,z), 
where the terms correspond the 1-halo (Poisson) and 2-halo (clustering) contribution. They are given by 



Pi(k) 



f°° dn 
/ dM -7T7 


p(fc,M)l 


/o dM 


P 



and 



dn , s p(k,M) 

dM -rv-MM)— 

dM y J p 



I 2 



where the radial Fourier transform of the density profile is given by 



p{k,M) = 4?r / dr r 2 p(r,M) 



sin kr 
kr 



(12) 



(13) 



(14) 



(15) 



Conveniently, p(0) = M. As noted by Seljak p9[|, the fact that on large scales P(k) must approach Pn n {k) imposes 
the non-trivial constraint 



- / dM —MUM) = \ 
pJo dM v > J 



dv f{v)b{v) = 1, 



(16) 



which formally holds in the case a(M, z) — > +00 for M —> 0. As for the mass function, we practically force this 
normalization to be numerically exact by adding a constant to b(M) in our smallest-mass bin. 
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The res ulting power spectra at different redshifts are shown on Figure [l] for the ACDM model (with parameters 
listed in § [ V B ) . Also shown are the power spectra from the Peacock & Dodds fitting formula. As noted by 
p2| , the agreement is good for all redshifts. Note that the only tuning involved is that of the functional form of the 
compactness parameter c(M,z) (Eq. (HJ). 




0.01 0.10 1.00 10.00 100.00 
k (h Mpc" 1 ) 

FIG. 1. Power spectrum of the DM density at different redshifts for the ACDM model. The solid, dashed and dotted curves 
correspond to the Halo model, the Peacock & Dodds fit, and the simulations, respectively. The thin horizontal dotted line 
corresponds to the Poisson noise in the simulations. 



B. Virial Temperature 



To model the baryons, we first consider the virial temperature of a halo which is given by 

where r„ iP hy S is the virial radius in physical coordinates and A c (z, f2 m , J7a) is the average virial overdensity which 
can be evaluated using the fitting formulae of |j3]. The value /j, = 0.59 (the number of particles per proton mass) 
corresponds to a hydrogen mass fraction of 76%. The factor j3 v is exactly equal to 1 for a truncated singular 
isothermal sphere. In the general case, it can be considered as an unknown normalization parameter, whose exact 
value is determined using numerical simulations. Values of (3 V between 1 and 1.2 provide good fits to numerical 
simulations p9|j40| . Here, we adopt /3„ = 1. 

Taking the gas temperature to be equal to the virial temperature, we can compute the density-weighted temperature 
(Eq. ||) of the gas which is given by 

_ l f M max d 

T P = = dM ——MT V , (18) 

P JMmix, dM 

where M m j n and M m a X: are mass limits introduced to model the finite mass resolution and box size of our numerical 
simulation (see § I V C|) . Note that, to ensure proper normalization of the total DM density (Eq. 0), this mass range 



is not used to compute the power spectrum of the DM density (Eq. |12|]). The resulting redshift evolution of this 
temperature is plotted on Figure |^, for different mass ranges. 
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C. Gas Profile 



To derive the profile of the gas pressure p g of halos, we further assume that the gas is in hydrostatic equilibrium 
and thus obeys 



dr 



Gp g (r)M(< r) 



(19) 



where p g is the gas density and M(< r) is the total mass within radius r. Approximating the total mass as that of 
the Dark Matter only, it is easy to show that, for the NFW profile (Eq. 



M(< r) = 4irp s rl 



ln(l 



(20) 



We need further assumptions to integrate Equation ( J19| ) for the pressure profile p g {r). As Cooray |21|, we consider 
the simplest assumption, namely that the gas is ideal (p g — j^~PgT g ) and isothermal within each halo, with a 
temperature equal to the virial temperature, i.e. T g (r) — T v . In this case, the hydrostatic equilibrium equation for 
the NFW mass profile yields 



with 



Pg( r ) = Pgov{r/r Sl A), v(x,X)=e A (l + x)-, 



A 



knT v 



(21) 



(22) 



This agrees with |2l| up to a different power for r s in the last expression. The normalization p g Q of the gas density 
profile is set by requiring that the baryon mass fraction is equal to that of the universe, i.e. by 



M g = Air 



dr r 



(r) = g-M. 



(23) 



Following pjj ], we avoid an abrupt cutoff at the virial radius by multiplying the profile function v{x) by an apodizing 
function which we choose to be 



E (^7^) ; where E(x) is the error function. 
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In the limit of small radii, the gas density approaches p g (r) oc e r » , where r g = ^ is the characteristic radius 
where the gas central density approaches a constant. It can be considered as a "core radius" in the gas distribution. 
The three characteristic radii r v , r s and r g are plotted as a function of mass M at z — and z = 3 on figure | The 
inner radius r g scales roughly as r s with mass, but tends to increase as z increases. This is due to the fact that, for a 
given mass, halos tend to have a higher temperature at high redshift (Eq. |l7]). In this model, the gas distribution is 
thus predicted to be more peaked at low redshift. 




10 J 
M (h M ) 

FIG. 3. Characteristic comoving radii for the halo model as a function of halo mass. The virial radius r v (solid line), the 
NFW scale radius r s (dashed lines) and the gas inner radius r 3 (dotted lines) are all shown at z — (thick lines) and z = 3 
(thin lines). Note that the curves for the virial radius are identical for all redshifts. It is the interplay between these radii 
which determine the behavior of the gas pressure power spectrum on intermediate scales. We indicate by arrows the comoving 
spatial resolution of the numerical simulation at each redshift. 



D. Pressure Power Spectrum 

As for the dark matter density (Eq. JT^]), the pressure power spectrum P p (k) is related to the pressure profile by 

P p {k,z) = P pl (k,z)+P p2 {k,z), (24) 
where 1-halo and 2-halo terms are given by 



P p i(k) = 



dM 



i\/„ 



dn 
dM 



p g (k,M) 



(25) 



and 



P P 2(k) = 



M ma: 
M mia 



dM -^ 6(M )M^) 



flin(fc), 



(26) 



where (M m j n , Af max ) is the mass range introduced in Equation (|18[). Here, p(k,M) is the radial Fourier transform 
(Eq. |l5|]) of the gas pressure profile p g (r). The mean gas pressure is p g — -^-p g T pi where the mean gas density is 

assumed to be p g = ^-p- The pressure power spectrum is plotted at different redshifts on Figure ||. 



G 
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0.01 0.10 1.00 10.00 100.00 
k (h Mpc" 1 ) 

FIG. 4. Power spectrum of the gas pressure at different redshifts. The thin smooth curves correspond to the analytical halo 
model (with Mmin = 5 x 10 10 /i _1 M© and A/ max = 8 x 10 14 /i _1 M©), while the thick broken lines correspond to the simulations. 



It is useful to consider the bias of the gas pressure relative to the dark matter density, 



It is shown for different redshifts on Figure ^. 




0.01 0.10 1.00 10.00 100.00 
k (h Mpc" 1 ) 

FIG. 5. Bias b p (k) of the gas pressure with respect to the DM density (see text) at different redshifts. As before, predictions 
for both the halo model (with M m i n = 5x 10 10 /!. _1 M© and Af max = 8 x 10 14 /i -1 M©) and for the simulations are shown. 



On large scales (k — ► 0), the 2-halo term dominates and P(k) ~ P\i n (k), so the pressure bias reduces to 

1 r M ™ x dn 

b p (k ^0)~— d M ——b(M)MT. (28) 

PT P J M^rn dM 
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In this limit, the pressure bias is thus simply the pressure-weighted average of the halo bias b(M). The large scale 
pressure bias (at k = 0.06ft. Mpc^ 1 derived from Eq (27]]) is plotted in Figure ^| for different mass ranges (M m i n , 
M max ). 
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FIG. 6. Large scale bias b p of the gas pressure as a function of redshift. The predictions from the halo model (for different 
mass ranges) and simulation are shown for k — 0.06/i _1 Mpc. 



IV. SIMULATIONS 



The halo model we described above relies on strong assumptions, namely that the mass density and pressure can be 
described statistically as a collection of spherically symmetric, isothermal halos in hydrostatic equilibrium. Although 
complicated structures, like walls and filaments, are known to exist in the hierarchical clustering picture, it remains to 
be established whether their influence on the resulting pressure power spectrum is sufficient to invalidate our simple 
halo model. In the following, we describe the high-resolution hydrodynamical simulations we have used to check the 
domain of validity of our analytical approach. 



A. Cosmological simulations using Adaptive Mesh Refinement 



Adaptive Mesh Refinement (AMR) is now widely used in astrophysics 41 - 44 48] and offers an interesting alternative 



to the Smooth Particle Hydrodynamics (SPH) algorithm. AMR was recently adapted to the cosmological framework 
by several groups |^2|,^9|. The main advantage of this grid-based method, compared to the well-known Smooth Particle 
Hydrodynamics (SPH) method pq j, is its ability to capture shock-waves and contact discontinuities using only a few 
resolution elements (typically 2-3 cells). In comparison, SPH demands at least 50 particles per resolution element, in 
order to minimize Poisson noise related to the discrete nature of the algorithm |47j]. 

The basic AMR idea is to refine recursively an initially coarse grid (defined here as the coarse level I = 0) with 
smaller and smaller cells of comoving size Ax — (£ > 0). These refinements take place in regions where interesting 
features of the flow demand better spatial resolution (see Figure ^). For example, in fluid dynamics, one usually 
refines sharp discontinuities like shocks or contact surfaces. In cosmology, refinements are used within high-density 
regions, in order to resolve properly the gravitational collapse of small scale halos and the progressive build-up of 
large, massive ones. Due to the highly compressible nature of the self-gravitating cosmological fluid, this property 
improves dramatically the quality of the numerical solution. 

In this paper, we use the results of a numerical simulation obtained using the newly developed RAMSES code, 
described in details in j23j. This code is based on general AMR techniques, with several specific properties that we 
briefly recall here. Contrary to patch-based AMR, where a small number of individual grids (say a few thousands) 
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define rectangular regions of higher and higher resolution (see for example |42|,[43|,[l8| ) , RAMSES is a cell-based AMR: 
each cell can be independently divided in sub-cells and thus can describe complex flow features of arbitrary geometry 
(see for example jl9],[32]]). The price to pay for this improved flexibility is a higher computational complexity: to 
describe the AMR grid, we use the "Fully Threaded Tree" data structure described by Khokhlov §52], where cells are 
placed in a recursive tree connecting each cell to its 8 sister cells and its 6 neighboring cells, with minimum memory 
overhead. 

The N-body module of RAMSES is similar to the Adaptive Refinement Tree (ART code) of Kratsov, Klypin & 
Khokhlov |Q, with some differences outlined in Teyssier p^| . The hydrodynamics module of RAMSES uses an 
unsplit, second-order Godunov scheme. Time integration proceeds in two different ways in RAMSES, either using 
a single, small time step for all levels of refinement or using nested, adaptive time steps in order to speed up the 
calculation for coarse levels. The last option was retained for cosmological simulations and turned out to speed up 
the calculation by a factor of 10, compared to the single-time step case. A complete presentation of the performances 
obtained by RAMSES is presented in Teyssier |23| , using standard tests to validate our algorithm and its application 
to cosmological hydrodynamics. 



B. Simulation Parameters 



In this paper, we present the results of a RAMSES simulation for a flat ACDM model with f2jv/ = 0.3, = 0.7, 
SIb = 0.039. We use the transfer function given in [ p4[ . The resulting power spectrum is normalized to the COBE 
data (5f| (this corresponds in our case to erg = 0.93y The comoving box size was set to 100/i _1 Mpc, imposing a 
maximum scale length in the simulation: k m i n = 0.06 h Mpc~ x . Initial conditions are specified using an initial grid 
of 128 3 dark matter particles. This corresponds to a mass resolution of M m ; n = 4 x 10 10 M©. The particles were 
initially displaced according to the Zel'dovich approximation up to the starting redshift Zi ~ 55. The baryon density 
and velocity fields were perturbed accordingly. We assume that baryons are described by a purely adiabatic, 7 = 5/3, 
fully ionized plasma. Since we are studying the SZ effect, we are mainly concerned by the pressure field of the baryons, 
which, at the scales of interest here, should remain unaffected (at least to first order) by neglecting others physical 
processes such as cooling, star formation and supernovae energy feedback. We plan to study the influence of these 
others physical inputs in a future paper. 

The key parameters of any AMR simulation are the refinement criteria used to dynamically create the refinement 
tree. We use the "quasi-Lagrangian" approach described in Kratsov, Klypin & Khokhlov E9|. The idea is to mark 
for refinement any cells whose gas or dark matter overdensity exceeds a level-dependent threshold given by 

^ = 0,5) = 1,80,640,5120,20480,163840 (29) 
P 

The coarse grid (£ = 0) is nothing but the 128 3 particle grid used to set up the initial conditions. The I = refinement 
criteria (p/p > 1) ensures that, initially, the whole computational volume is covered by a 256 3 grid {1 = 1). In this 
way, initial small scale perturbations are sampled by two grid points, which turns out to be necessary to avoid small 
scale power damping. At later times, large voids appear in the simulation where the resolution is locally degraded 
down to the 128 3 coarse grid. For higher levels of refinements, the density thresholds are chosen in order to refine 
cells that contains between 5-10 particles (or fluid mass elements). This "quasi-Lagrangian" approach preserves the 
initial resolution in physical coordinates, adapting the local comoving spatial resolution to collapsing mass elements 
(see j|5| for a complete discussion). 

The maximum level of refinement reached at the end of the present run is £ — 6, which gives a formal spatial 
resolution of 8192 3 (or 12 kpc /i" 1 comoving) at z = 0. The comoving spatial resolution scales roughly with redshift 
as (1 + We also use the adaptive time-stepping procedure of RAMSES to speed up the computation: 157 time 

steps only were necessary at the coarse level while 8185 time steps were necessary at the finest level. The Layzer-Irvine 
energy conservation was better than 1% at the end of the simulation. Starting with 19 x 10 6 cells at Zi = 55, we 
end up with 28 x 10 6 cells at z — in the tree structure (including the coarse level). Six different output times were 
considered: z — 0, 0.5, 1, 2, 3 and 5. 



C. Effect of finite mass resolution and box size 



The question that arises when interpreting numerical simulations is that of its range of validity. Here, we discuss 
three limitations of cosmological simulations which can introduce spurious effects: the finite box size, the mass 
resolution and the spatial resolution. 







Since we solve for the gravitational evolution of the system assuming periodic boundary conditions, the power on 
scales larger than the box size is completely suppressed. This translates into a maximum mass above which no halo 
can form. Indeed, the probability to find a large halo of mass, say, 1O 16 /i _1 M0 in our computational box is close to 
zero. An estimate for the maximum mass M max which corresponds to our box size is given by 

n(> M max )L box = 1 (30) 

where n(> M) is the cumulative Press-Schechter mass function. For Lbox = 100 /i _1 Mpc, this gives M max ~ 
7 x 1O 14 /i _1 M0. Note that, on large scales, we are dominated by the "cosmic variance" of the simulation: depending 
on the exact random field realization for the initial conditions, we may still form by chance a very large cluster in 
the simulation. The next step is to check the convergence of various quantities with respect to M max , using our 



analytical model. As we show below (§VE), we find that good convergence in the SZ power spectrum is achieved for 
M max = 2 x 10 1 /i _1 M0, or equivalently Lbox = 400 h~ 1 Mpc. On the other hand, for a fixed number of particles, 
increasing the box size inevitably degrades the mass and spatial resolutions of the simulation. Our current choice of 
Lbox = 100 /i _1 Mpc is a trade off between these two opposite constraints. 

The mass resolution, M m - m oc fc max , is related to the small scale power in the initial conditions. No object with 
M < M m i n will form, leaving instead a cold, smooth background between collapsing objects. In order to account for 
these two effects (finite box size and mass resolution) , we introduced in the analytical description a mass range where 
the different integrals are computed. This limited mass range is supposed to mimic the lack of large mass haloes due 
to the finite box size (M max ), and the lack of small mass halos due to the power cut off at small scale (M m - m ). 

The final numerical effect which needs to be examined is the spatial resolution. The spatial resolution is related to 
the minimum scale below which the code is not able to solve the equations. For example, in a standard Particle Mesh 
code, this scale is equal to 2 cell size of the underlying Cartesian mesh, below which the force between 2 interacting 
particles goes smoothly to zero. In practice, this acceleration cut off at small radii results in a damping of the power 
spectrum up to scales as large as 8 comoving cell size. For gas dynamics, the problem of spatial resolution is even 
more crucial. Shock waves need at least 2 cells to be properly resolved. If small scale velocity fluctuations are not 
properly sampled, no shock dissipation can occur at all. In this case, the small scale density fluctuations are only 
compressed adiabatically, i.e. they remain cold, even though their temperature should rise to the virial temperature. 
This can affect both the mean and the variance of the computed gas pressure. 

During the course of gravitational collapse, the initial perturbations contract at smaller and smaller comoving 
scales, making the spatial resolution issue even more crucial. The spatial resolution therefore translates directly into 
a minimum mass, below which halo collapse cannot be resolved properly by the code. For a Particle-Mesh code, 
coupled to a standard gas dynamics code, the resulting minimum mass can be as high as 100 particles, for a halo 



to be properly resolved. Most of the initial small scale power is then lost this way. As discussed in section § IV A , 
AMR techniques allows one to dynamically adapt the comoving resolution to the collapsing mass element. In light of 
the results obtained in this paper, we claim that the RAMSES code allows us to resolve properly small scale shock 
heating, down to the minimal halo mass imposed by the initial conditions. 



D. Computation of the Power spectrum 



In this paper, we need to compute the DM density power spectrum, together with that of the gas pressure. Since 
the spatial resolution goes from 96 kpc hT x at z = 5 down to 12 kpc hT x at z = for a box size of 100 Mpc h , the 
dynamical range imposes a severe challenge in computing power spectra. The brute force approach would consist in 
Fourier analyzing each field on a large Cartesian grid of 1024 3 for z = 5 (barely feasible on a supercomputer), up to 
a Cartesian grid of 8192 3 for z = (completely out of reach). 

A more appropriate approach is described in Jenkins et al. M] and Kravtsov & Kylpin |5Q ]. The idea is to define 
each field (dark matter density or gas pressure) in a nested structure of tractable mesh (we use 128 3 cells). At each 
level £, we divide the whole computational volume in £ 3 Cartesian 128 3 sub-cubes and add up the fields computed 
in all these sub-cubes together. We then use a standard Fourier power spectrum estimation for the resulting field. 
Each scale £ provides a reliable power spectrum estimation between 2 x [fc m ; n , fc max ], where fc m i n « 8(27r/Lbox) and 
&max ~ 16(27r/Lbox)- These limits are determined empirically, by comparing the results obtained for 3 nested 128 3 
grids and for a single 512 3 grid. The resulting power spectrum estimation is obtained by combining the different 
"band power" estimates together. Finally, we subtract from the dark matter density power spectrum the shot noise 
(or Poisson noise) due to the discrete nature of the particle distribution. 
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V. RESULTS 



A. Dark Matter Power Spectrum 



Figure [l] shows the dark matter power spectrum P(k) at different redshifts, derived from the Peacock & Dodds |3S[| 
formalism, the halo method and the simulations. The finite box size and the mass resolution (see the Poisson noise 
limit on the figure) restricts predictions of the simulations to the range 0.06 < k < 60 h Mpc -1 . For all redshifts, the 
simulations agree very well with Peacock & Dodds within this range, which spans 3 orders of magnitude. This range 
is consistent with the formal spatial resolution of the simulation, which spans 4 order of magnitude, after accounting 
for the small scale power damping discussed in §IVC. In comparison, the simulations used by Refregier et al. JToj ] 
and Seljak, Burwell, & Pen ||(| agreed with Peacock & Dodds only within one order of magnitude in k, and therefore 
seriously limited the SZ predictions of these authors 19 1. The present simulations therefore allow us to significantly 
improve over th ese pr evious works. 



As stated in § 111 A , the halo model prediction for P(k) is in good agreement with the Peacock & Dodds spectrum 
for all redshifts (see |32j). Small deviations of the order of a few percents can be noticed for k < 10 h Mpc -1 , but are 
acceptable for our purposes. Note in particular that the evolution of the break in the power spectrum between the 
linear and the nonlinear regime is well modeled by the halo model. 



B. Temperature Evolution 

Figure || shows the evolution of the mass- weighted temperature T p (z) (Eq. Q|) as predicted from the halo model 
(Eq. |ll|) and from the simulation. For the halo model, different mass cut off ranges M min < M < M max (see Eq. [ pL8| ) 
are examined. 

The halo prediction is sensitive to M max and M m [ n at low and high redshift, respectively. At high redshift, the 
characteristic non-linear mass scale M+ is close to the mass resolution of the simulation, for which no shock heating can 
occur. Logically, this translates in an underestimation of the mean temperature. At low redshift, the non-linear mass 
scale M+ approaches dangerously the maximum mass scale due to the finite box size. Here again, the temperature is 
underestimated, but this time, it is due to the absence of very high mass halos. 

A good overall agreement (mean temperature and power spectrum evolution) between the simulation and the 
analytical model is obtained for the following mass range: 

5 x 10 10 < M < 8 x 10 14 /i -1 M Q (31) 

It is reassuring that this best-match mass limits are close to the natural limits imposed by the simulation, as shown 



m 



IV C 



For a mass range of 5 x 10 10 < M < 4 x 10 14 h 1 M Q , the mean temperature evolution for the halo model (not 
shown on Figure S) is virtually indistinguis hable from the simulation curve, with T(z = 0) ~ 0.34 and 0.33 keV, 



respectively. However, as we will see below (§VE), the agreement with the pressure power spectrum slightly worsens. 
For M max = 8x 10 14 (our fiducial limit) and M max > 10 16 /i _1 M Q (the converged asymptotic limit) the halo model 
predicts T(z = 0) ~ 0.41 and 0.50 keV, respectively. 



C. Mean Comptonization Parameter 

Our calculation of the evolution of the density- weighted temperature allows us to compute the mean comptonization 
parameter (see Eq. Q|). For the simulation, we find y = 2.65 x 10 -6 . The halo model with our fiducial mass limits 
(Eq. pl|) predicts a value of y = 2.84 x 10 -6 . As expected, the agreement is slightly better for M max — 4 x 10 14 ft. _1 M Q 
with y = 2.72 x 10~ 6 . In all cases, these values are similar to those found in other studies (eg. jl4|,[l9 29|). It is 
interesting to study how our prediction for y depends on the mass limits. We find experimentally that 

/ M \ ~ - 06 / M \ °- 01 

( 8 xiffr-0 ■ < 32) 

Since length scales scale as Mi, this means that y oc A _0 - 2 _L 03 , where A is the effective resolution length and L 
is the box size of the simulation. This quantity is thus quite sensitive to the mass resolution, but not very sensitive 
to the box size. This should be kept in mind when interpreting and comparing the predictions of simulations. The 
asymptotic value corresponding, effectively, to an infinite mass range is y ~ 3.01 x 10~ 6 (for = 10, see ||). 



11 



D. Projected Map of the Pressure 



Figure [7| shows a map of the y-parameter projected through one side of the 100/i -1 Mpc box at z = 0. Clusters 
produces fluctuations in y of the order of 10 _4 -10 -6 , while supercluster filaments produce fluctuations of the order 
of 10 _6 -10 -8 . This figure is quite similar to the corresponding figure from Refregier et al. |ll| derived from the 
Moving Mesh Hydrodynamical code (MMH) of Pen |l^,|l8|. Note however, the more circular appearance of clusters 
in our simulations, compared to the more elongated cores seen in that of Refregier et al. [fl9|| . The number of small 
mass objects is also much higher in our simulation. The main difference between RAMSES and MMH is the spatial 
resolution. The improved dynamical range of RAMSES compared to MMH yields to the fragmentation of large scale 
filemcnts into small clumps. The effective mass resolution in RAMSES is consequently much smaller than the one 
reached with MMH. Moreover, substructures within large halos, clearly visible in Figure 0, are able to survive to tidal 
stripping much longer in RAMSES than in MMH. This explains the more elongated cores of large halos in MMH 
simulations. 



FIG. 7. Map of the y-parameter projected through on face of the 100ft 1 Mpc box at z = 0. 

Figure || shows a close up of one the clusters in the AMR simulation. Panel (a) shows the map of the y-parameter 
for the cluster, while the panel (b) shows the AMR grid for this region. Notice how the AMR grid has adapted to 
provide higher resolution for the cluster core and for the two substructures below the cluster. The virial radius for 
this cluster is about l-1.5/i _1 Mpc, and thus almost entirely fills the displayed region (3.12ft- 1 Mpc on a side). This 
figure also shows that the structure of individual clusters is more complicated than that assumed in our halo model. 
In the following, we show that the halo model nevertheless provides an acceptable prediction for the pressure power 
spectrum. 
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E. Pressure Power Spectrum 



Figure |I| shows the power spectrum of the gas pressure at different redshifts for both the fiducial halo model (with 
our fiducial mass range of Eq. |n]) an d for the simulations. The overall agreement is remarkable: the shape of the 
power spectra agree approximately for both methods, along with the scaling with redshifts which changes order above 
and below k ~ 3/i _1 Mpc -1 . According to both method, the pressure power spectrum does not evolve as much as 
the DM density power spectrum. Qualitatively, this is due to the fact that, at large redshift the smaller amplitude 
of the density power spectrum is compensated by the larger biasing of peaks with high temperature which dominate 
the pressure power spectrum. 

There are nevertheless quantitative differences: the halo model predictions have somewhat different normalization 
at low k and not exactly the same shape at large k. This can be seen more clearly by studying Figure ^, which shows 
the bias bp(k) of the gas pressure (Eq. |27| ]) at different redshifts. Again both the simulation and fiducial halo model 
are displayed. Note that we used the halo-model power spectrum in the denominator in both cases, so that Figure |^ 
is nothing but a change of scaling of the pressure power spectrum. 



FIG. 8. Close up on a cluster in the simulation, (a) Map of the {/-parameter for the cluster. The logarithmic color scale 
ranges from y — 10 -7 to 10~ 4 . The image size is 3.12/i~ 1 Mpc on a side, (b) AMR grid for this region with levels from 1 = 1 
to 6. 

The difference in shape at large k is more pronounced at low redshifts where the simulation pressure spectrum 
falls off faster than the halo model. The shape of the pressure power spectrum at intermediate and small scales is 
determined by the interplay between the different characteristic radii r v , r s and r g (see § 

at low redshifts the gas inner radius r g is smaller at the relevant mass (M» is about 10 14 /i _1 M Q ) and is dangerously 
close to the effective resolution limit (about 12/i _1 kpc) of the simulation. This would tend to reduce the pressure 
power at small scale in the simulation and thus explain the discrepancy. Numerical dissipation could also be at the 
origin of this small scale power damping (at a scale comparable to the spatial resolution). Indeed, this effect tends 
to induce a spurious increase in entropy, which results in a larger "core radius" in the gas distribution. Finally, the 
disagreement could also be caused by the limitations of the halo model. As Figure H illustrates, the structure of 
clusters is more complicated than that assumed in our model. In particular, cluster profiles may not be isothermal. 
In fact, our simulations indicate that the gas temperature tends to rise towards the center of clusters. This would 
lead to larger effective core radii and thus tend to suppress power on small scales. This effect could thus explain the 
discrepancy, as it is more pronounced at low redshifts. 

The mismatch of the normalization at small k can be studied using Figure [|, which shows the large-scale pressure 
bias as a function of redshift for both methods. The squares show the predictions for the simulations at the largest 
available scale, k = O.O6/1 Mpc -1 . The solid line shows the fiducial halo prediction at the same scale (derived from 
Eq. p7[). Both methods predict that the pressure bias increases with redshift, and agree relatively well for z < 1. 
The two methods however disagree by about 40% at z ~ 3. In addition, the change of slope predicted by the halo 
model at this redshift is not observed in the simulations. 

The large scale pressure bias depends on the mass range chosen for the halo model. This can be seen by examining 



III C). Figure shows that 
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the dotted curve in Figure It shows the halo prediction for a wide mass range, which has reached convergence to 
the asymptotic limit of an infinite mass range. This curve and that corresponding to the fiducial mass range agree at 
intermediate redshifts (1 < z % 2), but differ both at lower and higher redshifts: increasing M max tends to increase 
b p at low redshift, while decreasing M m ; n tends to decrease b p at large redshift. While this dependence must be taken 
into account, we could not find any reasonable mass range which eliminated the discrepancy with the simulations. 

The halo prediction for the asymptotic case of k — > (Eq. [p8[ ) differs from the k = 0.06ft. -1 Mpc curves shown in 
Figure [| only by about 20%. It therefore captures the dominating dependence of the large-scale pressure bias, and 
can thus be used to understand the origin of the discrepancy. For a given choice of the mass range, the asymptotic 
halo prediction is very robust. Indeed, once the mass function has been fixed to ensure an agreement with the DM 
density power spectrum, b p (0) only depends on the halo bias function b{v). After experimenting, we found that a 
modification of the Mo & White (37| relation for b{v) (Eq. |ll|) can indeed reduce the discrepancy. This is however 
poorly motivated, and would not be necessarily consistent with the observed clustering of DM haloes in N-body 
simulations p7pq -p9|]. The Sheth, Mo & Tormen Q b(v) relation was shown to better agree with N-body results 
|p9| ; its use in place of the Mo & White relation however slightly worsens the discrepancy for the gas pressure. We 
therefore do not pursue this possibility here. 

Instead, we note that this discrepancy for the large-sale pressure bias could come from the limitations of the halo 
model. It indeed assumes that that all the matter can be modeled as a collection of collapsed halo. At early times 
and on large scales, filaments and sheets indeed dominate the large-scale structure and are awkwardly described by 
the halo model. Visual inspection of SZ maps obtained by the simulation at large redshifts (z = 2-3) confirms this 
tendency qualitatively. 

Our results for the pressure bias are in good qualitative agreement with Figure 6 in Refregier et al. JlSj] . They 
however find a faster fall off at large k and a larger evolution of the pressure bias b v (k). These two facts can be 
explained by the coarser spatial resolution of their simulations, which, as discussed in §[VC, translates into a larger 
effective mass resolution (M m - in ~ 5 x 10 12 h~ 1 M & ). 



F. Sunyaev-Zel'dovich Power Spectrum 



From the estimates of the evolution of the pressure power spectrum and of the density- weighted temperature, we can 
compute the SZ power spectrum, using Equation (||) . The resulting spectra are shown in Figure ^ for the simulations, 
the halo model with our fiducial mass range, and the asymptotic halo model. The primordial CMB power spectrum 
for our ACDM model is also shown, as derived using CMBFAST [^0| . 

The agreement between our fiducial halo model (dashed line) and the simulations (dotted line) is excellent, showing 
that some of the discrepancies in the pressure power spectrum cancel out when integrated along the line of sight. 

The dot dashed line in Figure ^| shows the halo prediction for a wide mass range (corresponding to a good approx- 
imation to the asymptotic limit < M < oo). This model agrees with our fiducial model (and the simulations) for 
I > 7000, but exceeds it by a factor of as much as ~ 2 below that. Note that the difference comes primarily from the 
larger M max value, while a change in M m i n has very little effect. This shows that the main limitation of the simulation 
is the finite box size, which limits the number of massive clusters at low redshifts. 

The impact of spatial resolution can be studied by examining Figure |l0|. For this figure, the SZ power spectrum was 
computed from the fiducial halo model, after restricting the integration of the pressure power spectrum to different 
upper bounds A: max for the wavenumber k. (In all cases, the minimum bound was chosen to be fc m in = 0.06/i Mpc" 1 , 
the value imposed by the finite size of the simulation box). The figure shows that the SZ power spectrum is quite 
sensitive to the spatial resolution (i.e. to fc max ). A degradation of the spatial resolution to /c max ~ 5ft Mpc -1 leads 
to a severe drop of power for t > 2 x 10 3 (see also ]l9|). Interestingly, the SZ spectrum does not vary much if k max is 
taken be 20 ft Mpc -1 or above. For our simulation, /c max ~ 60ft. Mpc -1 at low redshifts, which is sufficient to reach 
convergence. The main limitation of our prediction from the simulation is thus that arising from the finite box size. 
On the other hand, thanks to the AMR scheme, spatial and mass resolution are not a limiting factor. 

Our prediction for the SZ power spectrum agrees approximately with the results of da Silva et al. up to a 
small overall rescaling due to our larger value of er 8 . Our results also agree with the prediction of Refregier et al. [^9| , 
within their £-range of confidence. As they noted, their predictions were limited by the spatial resolution of their 
simulations, which corresponds to /c max m 5ft. Mpc -1 . Our curve in Figure [w| corresponds to this value of fc m ax and 
is in good agreement with their result. 
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FIG. 9. SZ power spectrum for the halo model and for the simulations, in the RJ regime. The spectrum for the halo model is 
shown for different mass ranges: our fiducial mass range (dashed) and the asymptotic halo model (dot-dashed). The primordial 
CMB power spectrum for the A.CDM model is also shown for comparison. 
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FIG. 10. SZ power spectrum for the halo model with different values of fc max . Our fiducial mass limits of 
M min = 5 x lO lo /i _1 M and Af max = 8 x lO 14 ft _1 A/ , and a value of fc min = 0.06/i Mpc -1 were adopted. The power 
spectrum of the simulations is also shown for comparison. 



VI. CONCLUSIONS 



We have studied the 2-point statistics of the gas pressure and of the resulting SZ fluctuations using two methods: 
a RAMSES numerical simulation, using Adaptive Mesh Refinement, and an analytical halo model. Overall, the 
agreement is good, once the mass resolution and finite box size of the simulations are accounted for. The agreement 
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is almost surprising, given the simple nature of the halo model. The halo model plj is self-consistent and relies on 
the assumption of virial and hydrostatic equilibrium of the gas. After it has been tuned to produce the correct DM 
density power spectrum, the model does not require any new free parameters and thus has a strong predictive power. 

The predictions for the SZ power spectrum from our numerical simulations are mainly and solely limited by the 
finite box size. The AMR techniques indeed gives a reliable dynamic range of about 3 orders of magnitude for the 
power spectra, corresponding to a formal resolution of 8192 in linear scale, relieving us from the limiting effect of 
finite spatial resolution. The halo model allows us to extrapolate to an infinite box size, thus providing us with an 
accurate SZ power spectrum in the full multipole range 10 2 < t < 10 . This is an important improvement upon 
earlier calculations [p9|j20| which were severely affected for I > 2 x 10 3 by their coarser spatial resolu tion. 



In principle, the effect of finite box size could be alleviated using larger simulations, as discussed in §[V_C. However 



to achieve convergence, one would need to reach a limiting mass of M max = 2 x 10 15 ft. _1 M Q /i _1 which would correspond 
to a box size of 400 Mpc in a ACDM model. Given our current computer resources, this is prohibitive for the 
near future. One thus has to rely, as we have done, on analytical models to extrapolate beyond the sampling variance 
of the simulations. Analytical models, like the one we have presented here, also have the advantage of providing a fast 
way to explore a wide range of cosmological parameters. They also help to provide a physical understanding of the 
physics and scales which contribute to the SZ statistics. Using our halo model, we have shown for example that two 
characteristics length scales contributes to the SZ power spectrum: the virial radius r v and the gas inner radius r g . 
It would be interesting to understand in more details how the SZ power spectrum provides a measure of the different 
regions of the NFW profile on a range of mass scales. Another interesting avenue is the study of processes such as 
feedback and cooling on the SZ statistics. This will be explored in future work. 
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